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Abstract 

The self consistent version of the density functional theory (DFT) is presented, which 
allows to calculate the ground state and dynamic properties of finite multi-electron systems 
such as atoms, molecules and clusters. The exact functional equation for the effective 
interaction, using which one can construct the action functional, density functional, the 
response functions, and excitation spectra of the considered systems, are discussed. We 
have also related the eigenvalues of the single-particle Kohn-Sham equations to the real 
single-particle spectra. We begin with the standard action functional and show that it is 
useful in calculating the linear response functions \ via which all physical characteristics 
of a many-body system can be expressed. It is shown, that function x can De causal 
(retarded), noncausal and advanced ones. This resolves the well known paradox related 
to the causality and symmetry properties of the response functions and the effective 
interaction. 
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INTRODUCTION 

The density functional theory (DFT) originated from the pioneer work of Hohenberg and Kolin 
about thirty five years ago. Since then an impressive progress has been achieved by using 
DFT in describing microscopic many-particle systems and especially in studying multi-electron 
objects, such as atoms, molecules clusters and fullerenes. The recognition of DFT successes 
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culminated in awarding W. Kohn the Nobel Price in Chemistry in 1998. At the beginning 
DFT was limited to consideration of only the ground state properties of many-particle systems, 
leaving aside their dynamic properties, which are closely related to the system's behavior in the 
time-dependent external fields. This limitation of DFT was overcome by Runge and Gross [0] 
who thus have transformed DFT into the so-called time-dependent density functional theory 
(TDDFT). Both, DFT and TDDFT are based on the one-to-one correspondence between the 
particle densities of the considered systems and external potentials acting upon these particles. 
For the sake of simplicity and definiteness, we assume that the time dependent part of the 
considered electronic system's density p(r,t), which is created under the action of an external 
time-dependent field Xv ext (r, t), is developing from the system's ground state. As consequence 
of it one can conclude that each observable, describing the system, can be written as a functional 
of the density, among them are the functionals of the ground state energy E[p] and of the action 
A[p] |l|, 0. Then DFT and TDDFT establish an exact correspondence between an interacting 
many-body system and a fictitious non-interacting Kohn-Sham system [JTJ, £| E|. The main 
advantage of the mapping is that Kohn-Sham system is described by a system of Hartree- 
type single particle equations, which is not too difficult to solve. As a result, after solving 
these equations in the case of an external time-independent potential, one can obtain and 
predict, at least in principle, the atomic, molecular, cluster and solid bodies binding energies, 
phonon spectra, activation barriers, and etc (see e.g. jT]). The same is true in the case of 
a time-dependent potential when the solution, generally speaking, yields the single-particle 
and collective excitation spectra, describes the behavior of a system in strong external fields, 
including that of atoms in short laser pulses (see e.g. [||, [|]). Unfortunately, the one-to-one 
correspondence establishes only the existence of the functionals in principle, leaving aside a very 
important question: how one can construct them in reality. This is why the success of DFT and 
TDDFT strongly depends upon the availability of good approximations for the functionals. As 
such, an approach based on the local density approximation and on numerical simulations of 
the ground state energy of a homogeneous in density electron system was suggested already in 
However the real molecules or atoms are not homogeneous and numerical simulations are 
not universally good. So, both justification and improvement of this approximation is required. 
But attempts to do it run into hard problems. They had to be overcome and it is highly 
desirable to find a systematic way to construct the required functionals. On the other hand, an 
important element of DFT, which is the action functional A[p], is suspected to be ill-defined 
because of a number of contradictions in it. Perhaps the most important among them is the 
contradiction, which came from the analysis of the causality in TDDFT. It appeared, that this 
contradiction play a central role in creating a number of other difficulties |], [6[]. In fact, it 
signals the existence of a very serious paradox: on one hand, according to the Runge-Gross 
theorem |2[] A[p] is a well defined functional, while on the other hand, the use of this functional 
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leads to fundamental difficulties. 

Here we report some recent developments in eliminating this paradox. As we shall see, 
the methods used in resolving it give a good opportunity to construct equations defining the 
functionals A[p], E[p] and related to them effective interaction and the linear response function. 
Thus, this Comment will focus on formulating such an approach to the self-consistent density 
functional theory, that is based not only on the theorems of the existence of the functionals but 
also on the exact equations for them. 

We will start with discussion of the paradox and show how it is eliminated, then go to 
single particle potentials, present the functional equations for the energy and action, derive 
the expression for the effective inter- electron interaction and illustrate the effectiveness of the 
suggested scheme. 

It will be demonstrated that the presented approach is at the same time rather simple, as 
compared to the usual, say, diagram many-body technique, and offers a possibility to calculate 
relatively easy all the required characteristics of any multi-electron systems. 

SYMMETRY AND CAUSALITY 

Let us start with the discussion of the action functional A[p], which is of key importance for 
the subject of this paper, because it permits to express the reaction of the many-body system 
under consideration to an external field. Due to the existence of one-to-one correspondence 
between the external potential Xv ext (r, t) [Q] and the time- dependent density p(r, t), the external 
potential determines also the time-dependent wave function 4>[p](t) of the system in question. 
Therefore, the expectation value of any quantum mechanical operator is a unique functional of 
the density . This statement is also true for the action functional A[p], which is determined by 
the expression [|, |J, 



The values of ti and tf can be chosen arbitrarily, we take t« = —oo;tf = +00. The time- 
dependent wave function (j)(t), which corresponds to a stationary point of the action functional 
(1), is a solution of the Schrodinger equation, 




(1) 



i-<j>{t) = H(t)<f>(t) 



with the hamiltonian, that we present as a sum of two terms, 



H(t)=H + H l (t). 



(2) 
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Here H = K + V acting alone describes the ground state of the considered system, while Hi(t) 
represents the external field and is given by the following expression: 



H 1 (t) = J p(r)Xv ext (r,t)dr. 



Here p(r) is the electron density operator. The operator K is the hamiltonian of non-interacting 
particles and V represents the two-particle Coulomb interaction. Upon using eqs. (1,2) the 
functional A[p] can be presented as, 

A[p] = A 1 [p] - J Xv ext (r,t)p(r,t)drdt, (3) 

with A\ = A\\=o being a universal functional of the considered system, which is explicitly 
independent upon the external field. Using the action functional A\, it is possible to construct 
the inverse of the linear response function Xi i- e - X" 1 ? [U> 

X~ l r l5 1 1? r 2 , t 2 ) = 1[H 4 

5p{r 1 ,t 1 )5p{r 2 ,t 2 ) 

One needs to know the function x, because all characteristics of a many-particle system are 
expressed via it. From eq. (4) it might seem as self evident, that not only x~ l but also the direct 
linear response function x defined by x _1 ( r i) ti, r 2 , t 2 ) must be a symmetrical function of its 
space and time arguments |3], |5|, |J. However, x, being symmetric under interchange of (i"i,ti) 
and (r2,t2)j cannot be a causal function, while determined in the usual way x( r ij r 2, ^2) is 
known to be a retarded (causal) function and, therefore, must vanish t 2 > t\ 0, ||. Therefore 
it was concluded that eq. (4) is incorrect ^ because the symmetry and causality looked 
as contradicting each other. To overcome this difficulty, it was suggested that the problem is 
related to the boundary conditions 5<p{tj) = 5<p(tf) = which one has to enforce on the variation 
of the wave function while deriving the time-dependent Schrodinger equation from eq. (1). It 
was argued, that since the action functional is defined only on the sufficiently narrow domain of 
v -represent able (VR) wave functions and densities, the condition 5(j)(tf) = cannot be fulfilled. 
The v -represent ability means that the wave functions and corresponding densities are solutions 
of the Schrodinger equation with an external potential v(r, t) [0, Thus, a non- vanishing 
boundary term appears when the variation 5A is made. Therefore, it was recommended to 
search for another functional to avoid the mentioned above causality-symmetry problems 0. 
But, as it was demonstrated above, the existence of the action functional is guaranteed by the 
Runge-Gross theorem 0], which play the key role in establishing TDDFT. As a result, the 
paradox has thrown doubt on the reliability of TDDFT and it became crucially important to 
resolve it. 
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In eliminating this paradox, we shall follow the papers M, |10[. Consider a stable multi- 
electron system evolving from its ground state 0o with the density po under the influence of a 
weak external potential Xv ext , that is in the limit A — > 0. This condition prevents the system 
under consideration from a steady heating, ensuring the existence of the linear response function 
X as well as its inverse x -1 [0? HI- F° r suc h a system all the densities and wave functions in 
the vicinity of the ground state are VR including those which satisfy the boundary conditions 
H 11, |12|]. Therefore the correct time-dependent density can be obtained from the Euler 



equation 0, 

W'.t) 

It is worth to note, that this equation is more general than the boundary condition. It means 
that (5) can be correct even if S(fi(tf) ^ 0. Now we can use (5) to prove (4). According to (3), 
the Euler equation (5) can be written in the following form, 

SA[ P ] _ 



5p(r,t) Sp(r,t) 

To prove (4), we expand the solution of (6) in terms of the density fluctuations 5p = p\ — po, 
induced by the external field: 

6Al \ l 52Ai \ x \ \ n (7\ 

~Jp\ p0 7p2~' Po P _ Av «etU-* = °> (') 

where p\ is the density in the presence of Xv ext . The first term in (7) vanishes because of 
the stationary condition (5). Upon using the definition of the linear response function \ = 
Sp/S\v ext , the exact equation is obtained [|J, 

x -\nM^t 2 ) = , , 6 , A ' [p j , J a= . (8) 

dp(r 1 ,t 1 )dp(r 2 ,t2) 

Thus, we come to the conclusion that the second functional derivative of the action integral Ai 
appeared in eq. (8) defines and is defined by the inverse of the linear response function, which, 
as it should be, is a symmetrical function of its variables. To understand the physical meaning of 
eq. (8), which defines a symmetrical function consider eq. (1) which determines the action 
A. It is known that by varying this equation with respect to <p* one derives the Schrodinger 
equation for the retarded solution 0, <fi r (r,t), while the variation with respect to 4> leads to the 
advanced solution 0*, 0*(r, —t). As a result, the variation of the action functional with respect 
to the density p produces a superposition of the advanced and retarded solutions. This is why 
eq. (8) defines a symmetrical function. Now we are going to demonstrate that having the 
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noncausal function x _1 , it is possible to construct any type of the linear response functions x : 
noncausal, and causal (both retarded and advanced) linear response functions. The response 
function x is connected to its inverse by the following relation: 

5{t x -t 3 )<J(r! - r 3 ) = / x (r 1 ,v 2 ,t 1 ,t 2 )x(r 2 ,r 3 ,t 2 ,t 3 )dr 2 dt 2 , (9) 

We see that noncausal x" 1 is the kernel of the integral equation (9). The character of solutions 
X obtained from (9) depends on the boundary conditions which are defined by numbers T\ and 
r 2 . If the linear response function has to be causal, that is to be nonzero only at t 2 > t 3 then 
n = £3 and t 2 = 00. If noncausal x is required, then n = — 00 and t 2 = 00, while in the case of 
advanced x it has to be t\ = —00, and r 2 = t 3 . As we will see below, using a different way, the 
linear response function may be chosen indeed causal, noncausal or advanced. And vice versa, 
starting, for instance, with the causal linear response function x, taken as the kernel of (9), 
one can construct x~ l of desirable kind using the same eq. (9) and choosing the appropriate 
boundary conditions. Thus, we come to the conclusion: it is the existence of the symmetrical 
(noncausal) inverse x _1 that leaves no room for the paradox 0, |K|. Note that the particular 
case when x _1 is noncausal while x is causal has been studied [13| recently. 



It is worth to consider in more details the case when both x an d X -1 are causal since one 
may suspect that there are no such solutions of (9) [13|]. To show that the causal function x^ 1 
does exist we consider the well known definition of the linear response function, 

rti 

5p(r 1 ,t 1 ) = x( r i,r 2 ,t 1 ,t 2 )Xq(r 2 )v(t 2 )dr 2 dt 2 . (10) 
J 

For simplicity we substitute the external field in eq. (10) by v ext (r, t) = Xq(r)v(t). Then, eq. 
(10) transforms into Volterra's integral equation of the first kind, which always has a solution, 
see e.g. [fL4|. Thus, we conclude that the causal inverse Xr 1 exists. 



In order to demonstrate that one can use the causal (retarded), advanced and noncausal 
linear response functions equally effective let us consider the second functional derivative of 
the action with respect to external field. Because of the one-to-one correspondence between 
external potential Xv ext (r, t) and time dependent density p(r, t) such a functional Ai[v] exists. 
Then, taking into account that the inverse of the linear response x is given by 5Xv ex t/Sp = x" 1 
and using (8), one obtains, 

^Ailv] . . . 

X(ri,ti,r 2) t 2 ) = tt - ( tttt 7 — ttU=o- (11) 

dXv ext {r 1 ,t 1 )dXv ext {r 2 ,t 2 ) 

It is useful to derive eq. (11) directly from the definition (1). We restrict ourselves to an 
infinitesimally weak time-dependent field putting A = after completing the calculations. 
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Therefore the usual many-body technique of diagrams can be applied in calculating the varia- 



tional derivatives with respect to the external field []T5], [16|]. The variation AAi of the action 
Ai[v] induced by the time-dependent external field is then defined as, 



AAi 



dt{</>(t) 



4-H 
dt 



*f <0o|T(g rl (t)^(oo))|0 o > 
<0o|T(S(oo))|0 o > 



;i2) 



where T is the usual time-ordering operator, and H T i(t) is Hi(t) in the interaction representation 
with respect to the hamiltonian H. The expression for S(oo) is defined by the following way, 



S{po) = Texp 



H Tl (t)dt 



(13) 



where c is the contour of integration, which starts at — oo and ends at +00. Note that the 
functional Ai[v] formally defined as, 



A x [v] = iln< <p \T(S(oa))\(f) >, 



(14) 



will give the variation AAi with respect to the external field, which coincides with the variation 
following from eq. (14). Calculating the second variational derivative of Ai given by eq. (14) 
with respect to the external field Xv ext and putting A = in the final result, one gets the 
following relation: 

$ 2 Ai[v] . . 

11 J 1 -x(ri,*i,r 2 ,* 2 ) (15) 



5Xv ext (ri, ti)5\v ext (r 2 , t 2 ) 
-i(< </) I Tp(ri,ti)p(r 2 ,t 2 ) | <f>o > - < <t>> 



a=o- 



P(ri,*i 



&o >< 



5»o > 



It is seen from eq. (15) that the second functional derivative of the action integral A x with 
respect to the external field is defined by the noncausal linear response function, which, as it 
should be, is a symmetrical function of its variables [ffl. After performing the Fourier transfor- 
mation with respect to time, x( r i5 r 2,^) takes the form, 



X 



E 



< I pt(ri) I k >< k I p(r 2 ) | > < | p t (r 2 ) | k >< k \ p(n) | > 
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(E k - E ) + i5 



u + (E k - E ) - iS 



(16) 



Here | >, E and | n >, E k are the many-particle system's exact ground and excited state 
wave functions and energies, respectively. It is evident from (16) that \ has poles in // and 
IV quadrants of the complex u plane. It contains exactly the same information related to the 
properties of a system as the causal linear response function, while both these functions coincide 
on the imaginary axis of the plane, being related one to another by analytical continuation 
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H [RJ. Just in the same way the response functions of higher orders can be constructed. So, 
having at hand the action one can study the ground and excited states of a system. Thus, 

eq. (8) is of key importance in calculations of the excited states within the TDDFT framework 



Note, that eq. (8) was later verified in [20]. It is pertinent to point here that one can use the 



non-equilibrium Green function theory |17|, [18| to construct both the retarded (causal) linear 
response function \ri which has poles in the III and IV quadrants, and the advanced linear 
response function Xa, with poles being in the / and II quadrants. The only formal difference 
from the used above technique is the appearance of the counter ordered evolution operator T cr 
instead of the usual time-ordering operator T Thus we can use eq. (14) to evaluate A\[v] 



replacing T by T cr , in this case eq. (14) is modified as, 

Ax[v\ = iln < (j) \T cr (S(oo))\(f) >, (17) 

with T cr being the contour-ordering operator, which orders the operators along the contour cr 
according to the position on the counter of their time arguments. Then, upon using eq. (17), 
we have instead of eq. (15), 

5 2 A 1 \v] 

7T 2 |a=o= Xr(ri, h, r 2 , t 2 ) = -i0(*i - t 2 ) < 0o I [p(ri, h)p(r 2 , t 2 )] \ 4>o > ■ (18) 

Here the brackets [...] denote the commutator of the density operators. In the same way as 
Xr was constructed, one can get the advanced linear response function Xa, choosing a special 
contour ca along which the contour-ordering operator T ca orders the operators. Thus, different 
forms of linear response functions correspond simply to different choices of the contours. This 
result is quite obvious: if one marks the moment on the contour then it is possible to construct 
the retarded (advanced) linear response function, using eq. (17) for the action functional. If 
no moment of time is marked along the contour, and we deal with the original action func- 
tional given by eq. (14), we shall have the linear response function x defined by eq. (15) 
in accordance with the previous discussion. Thus, we can conclude that functional A is in 
fact well defined, while eq. (4), which determines the inverse function opens the possi- 
bility to calculate the dynamic properties of a system including its excitation spectra 0, IT3, [Z0 |. 



SINGLE PARTICLE POTENTIAL 

Now let us construct the single particle potential v s [p](r,t) of the time- dependent Kohn-Sham 
equations. It is convenient to define the exchange-correlation functional A xc 0. [| . 

A xc [p] = A [p] - AM ^^^ dtdrxdv,, (19) 

Z J I"i — Y 2 \ 
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in atomic units, used throughout this paper. Here A® is the functional of non-interacting 
Kohn-Sham particles, 

A [ P ]= f tf dtmt)[ P ] -° 

Hi 



*(*)[/>]}, 



with $(t)[p] being the unique time-dependent Slater determinant of the density p 0, |. The 
effective interaction is defined as f|, [19], PR , 

R^(r utl ^t 2 ) = + ff- (20) 

|r x -r 2 | dp(Ti,ti)dp(T2,t 2 ) 

Combining eqs. (8) and (20), the effective interaction R(°>(ri, t±, r 2 , t 2 ) is determined by the 
relation, 

i*(°>(n, ti, r 2 , £ 2 ) = xo '(ri, ti, r 2 , t 2 ) - x" 1 ^, h, r 2 , t 2 ), (21) 

where Xo 1 — <5 2 A)/<5p 2 is the inverse of the linear response function of non-interacting Kohn- 
Sham particles. In the same way, calculating the functional derivatives of higher orders, one 
can obtain the effective interaction R"\ I > 2, 



oil) 1 4- 4- 4- \ <^' +2y 4xc[p] 

ti (ri, ti, r 2 , t 2 , ...r/ +2 , ti +2 ) — —, 7~T7 / -r — 7 — \ 

bp{r 1 ,t 1 )dp{v 2 , t 2 )...5p(r l+2 , t l+2 ) 



A=0 



(22) 



Then, taking the Fourier transform of the effective interaction with respect to time, making 
the analytical continuation to obtain the causality, and upon carrying out the inverse Fourier 
transform, one can get the causal effective interaction terms 3?^. These terms IRy 1 are members 
of a series, which define the causal effective interaction 



M[p](r 1 ,t 1 ,r 2 ,t 2 ) (23) 

f 1 

= / 7[^ (/) ( r i^i ; r 2^2, ...T l+2 ,ti +2 )5p(T 3 , t 3 )...5p(r l+2 , ti +2 )dr 3 dt 3 . ..dxi +2 dt w . 

i>o L 

Having at hand one can construct the single-particle potential v s [p] (r, t) by integrating the 
equation, which determines this potential 

Sv[ P )(y) = Sv M (r tl) + 1 = 

<5p(r 2 ,i 2 ) 5p(r 2 ,t 2 ) |ri-r a | 

with f xc being the exchange-correlation potential. Eq. (24) can be integrated as simple as one 
can integrate the Tailor expansion of a function, 

VslpKruh) = J »(°)(r 1 ,t 1 ,r 2 ,t 2 )(J/>(r 2 ,t 2 )tZr 2 tZt a (25) 
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f 1 

+ H / 77— jry^ {0 ( r i^i) r 2, *2, -ri+2, £z+ 2 )<fy( r 2, t 2 )...5p(r l+2 , ti +2 )dr 2 dt 2 ...dr i+2 dti +2 . 



Then, one can use the exact functional equation for R(r l ,r 2 ,uJ, g) ||, [H| to construct approxi- 
mations to v s [p], 

R{r 1 , r 2 , u, g) = - — - — , (26) 
|ri - r 2 | 

1 ^ f f 9 ( , , ■ I A > dw A ' 

~oTl u~7 W / X{r 1 ,r 2 ^w,g)— -dri dr 2 — dg . 

2 op(ri, u)op{r 2 , —to) J Jo |rj_ — r 2 | 27r 

Here i2(ri, r 2 , cu, (?) is the effective interaction depending on the coupling constant g of the 
Coulomb interaction. The coupling constant g in eq. (26) is varied in the range (0 — 1). When 
g = 1 one has -R*- -* = R(g = 1). In eq. (26) the linear response function x( r i^ r 2,^, g) is the 
Fourier-image in (ti — t 2 ) of the linear response function. It is obtained from eq. (21) through 
algebraic transformation and by introducing the coupling constant g, 

X(*i, r 2,w,g) = Xo(ri,r 2 ,o;) + J Xo(r 1 ,r' 1 ,u)R(r' 1 ,r' 2 ,uj,g)x(r' 2 ,r 2 ,uj,g)dr' 1 dr' 2 , (27) 

with xo being the linear response function of non-interacting Kohn-Sham particles, moving in 



the single-particle time- independent field 0, [19]. It is evident that the linear response function 
X(g = 1) g° es to X as g goes to 1. 

The physical meaning of (25) is quite transparent: as it is seen from (24,25), the functions 
are directly obtained by calculating the corresponding functional derivatives of v s [p] with 
respect to the density. This shows that the single particle time-dependent potential contains 
a lot of information about the system. And, vise versa, to construct v s [p] one needs the same 
information. It is important to note that v s [p], determined by eq. (25), is suitable to describe 
the system in a strong external field due to the general property, given by eq. (3), with the 
provision that the series in eq. (25) is convergent. 

It is instructive to consider eqs. (8,15) in the case when external field Xv ex t is not weak. 
Deriving (8,15), we used the fact that the linear response function and its inverse exist when 
A = 0. In the case when X = X Q ^ 0, the existence of both of these functions is not self-evident. 
Moreover, it is quite probable that in the general case their existence is impossible to prove. 
The behavior of a system in an arbitrary (including the strength) time-dependent field may be 
very complicated and unstable, for instance, due to a permanent transfer of the energy from the 
external field to the system under consideration. On the other hand, the non-equilibrium the- 
ory is fitted to describe time-dependent steady states, it is those states that are stable against 



small perturbations [17]. But stability against small perturbations means that there exist the 



response function of the system under consideration and its inverse function. So let us suppose 
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the existence of both of these functions. This can be correct provided the series in eq. (25) 
are convergent. Under this condition we can apply the same consideration that has been used 
when deriving eqs. (8,15): the density of the system imbedded in the external field is VR, 
then the densities in the neighborhood are also VR fLlf . As the result, the action functional 
A[p] is defined, and eq. (5) is valid. One can also apply eq. (7) to prove (8), having in mind 
that A — > A , while p in the considered case is the time-dependent density in the external 
field XoV ext . On the other hand, we can use eq. (17) evaluating Ai(X ), which is correct in the 
presence of an external field. The symbol (Ao) emphasizes that the external field is not a weak 
one. As a consequence, we obtain the retarded response function Xr(Ao), given by eq. (17). The 
advanced response function Xa(-^o) is defined in the same manner by choosing the contour ca. 
The contour, which starts at — oo and ends at +oo, will produce the response function x(Ao). 
Thus, we are coming to the conclusion that the standard definition of the action functional 
(1), as well as eqs. (8,15,17,18) are also valid in the considered case of a strong external field. 
Further investigation of the behavior of a system in this case will be published elsewhere. 



EQUATIONS 



In this section we briefly outline the derivation of the functional equation for the exchange- 
correlation functionals E xc and A xc when system in question does not perturbed by an external 
field. In that case one has 

E xc [p] = A xc [p]\ p ( r)U1= o), (28) 

since A xc is also defined in the static densities domain. The exchange-correlation functional 
may be obtained from 0, 19|: 



A xc = -\ [ [x(ri, r 2 , iw) + 2>Kp(r 1 )5(w)S(r 1 - r 2 )] , — ~ -^dx 1 dT 2 - (29) 

2 J \Ti — r 2 | g 27r 

Eq. (29) presents the well-known expression for the exchange-correlation energy of a system, 
see e.g. 0, [7|, [19| . The only thing we need, in order to consider eq. (29) as describing A xc , is an 
ability to calculate the functional derivatives of A xc with respect to the density. This, as it is 
seen from eq. (29), reduces to the ability of calculating the functional derivatives of the linear 
response function \ with respect to the density p(r,u) which was developed in |jl9| , 22fl . Then 
the linear response function is given by (27), and the effective interaction R is defined by the 
exact functional (26). The single-particle potential v xc , being time-independent, is of the form 

v xc (r) = £p( r £\A a:c \p = p . (30) 
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Here the functional derivative is calculated at p = po with po being the equilibrium density. By 
substituting (29) into (30), it can be shown that the single particle potential v xc has the proper 

— 1/r, (TP], j23|. The potential v xc determines the energies 



asymptotic behavior, v 
Ei and wave functions < 



oo 





of fictitious Kohn-Sham particles, 



r H (r) + V ext (r) + v xc (r 



that in turn form the linear response function \0i 
Xo(ri,r 2 ,cj) = Yj 71 ^ 1 ~ w jO#( r i)&( r 2)0fc( r 2)0*( r i 



i.k 



and the real density of the system p, 



(31) 



p(r) = Y, n i\H r ) 



1 



1 



uj - uj ik + ir) lo + io ik - irj 



(32) 



(33) 



Here are the occupation numbers, V ex t contains all external single-particle potentials of the 
system, say the Coulomb potentials of the nuclei. Then, Vh is the Hartree potential, 



P(ri 



rdr h 



Vh(v) 



and ujik is the one-particle excitation energy uj^ = Si — Ek, while r\ being an infinitesimally 
small positive number. The described above equations (29-31) and (26,27) solve the problem of 
constructing the self-consistent DFT: one can calculate A xc , the ground state energy and exci- 
tation spectra of a system without resorting to approximations for A xc based on additional and 
foreign to the considered problem calculations such as Monte Carlo simulations, or something 
of this kind. We note, that using these approximations, one faces difficulties in constructing 
the effective interaction of finite radius and the linear response functions p, |5|]. On the base on 
the suggested approach one can solve these problems. For instance, in the case of a homoge- 
neous electron liquid it is possible to determined analytically an efficient approximate expression 
Rrpae for the effective interaction R, which essentially improves the well-known Random Phase 
Approximation jTJ| [T| by taking into account the exchange of electrons properly, thus forming 
the Random Phase Approximation with Exchange. The corresponding expression for Rrpae 
is as follows 

Air g 

R RP AE(q,g,p) = — + R E {q,g,p), (34) 
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where 

RE{q,g,p) = 

Pf 

Here the electron density p is connected to the Fermi momentum by the ordinary relation 
P = Pf/37t 2 . Having at hand the effective interaction Rrpae(q, g, p), one can calculate the 
correlation energy e c per electron of the electron gas with the density r s . The dimensionless 
parameter r s = r /a,B is usually introduced to characterize the density, with r being the average 
distance between electrons, and ap is the Bohr radius. The density is high, when r s < 1. 



12pl 



In 



1 - 



Ap 2 F 



2pf 
3q 



In 



2p F - q 



2p F + q 



1 

+ 3 



(35) 



Table 1. Correlation energy per electron in eV of an electron gas of density 

r s . The Monte Carlo results |24[| e* f are compared with the RPA 
calculations and with the results of [pi]] . Srpa denotes the results of 



RPA calculations, and Erpae denotes the results of the calculations 



when the effective interaction R was approximated by Rrpae [21 . 



r s 


e M 

c c 


£rpa 


£rpae 


1 


-1.62 


-2.14 


-1.62 


3 


-1.01 


-1.44 


-1.02 


5 


-0.77 


-1.16 


-0.80 


10 


-0.51 


-0.84 


-0.56 


20 


-0.31 


-0.58 


-0.38 


50 


-0.16 


-0.35 


-0.22 



As can be seen from Table 1, the effective interaction RppAE^q, p) permits to describe the 
electron gas correlation energy e c in an extremely broad interval of density variation. Note, 
that even at r s = 10 the mistake is no more than 10% as compared to Monte Carlo calculations, 



while the result becomes almost exact at r s — 1 and is absolutely exact when r s — > [pi . 

Now let us calculate the single particle energies ej, that, generally speaking, do not coincide 
with the eigenvalues £j of eq. (31). We remark that the eigenvalues do not make a physical 
sense, see e.g. To calculate the single particle energies one can use the Landau equation 

B, 

SE , . 

ST"- < 36 > 

Eq. (36) can be used since, as it follows from eqs. (27,32,33), the density and the linear 
response function depend upon the occupation numbers. Thus, one can consider the ground 
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state energy as a functional of the density and the occupation numbers, p5 



E[p(r), ni ] = T k \p(r),m] + ±J V H (r)p(r)d(r) + J V ext (r)p(r)d(r) + E xc [p{v),n t }. (37) 

Here T& is the kinetic-energy functional of non-interacting Kohn-Sham particles. It follows 
from eqs. (29,31,36,37) that the single-particle energies e« can be presented by the following 
expression, 



e~ < (j>i\v a 



> 



1_5_ 

2 5rij 



X(ri, r 2 , iw) + 2irp(r 1 )5(w)5(v 1 - r 2 ) 
ki - r 2 | 



dwdg'dridr 2 
2~n ' 



Here, as it follows from (27), 5x/8rii is given by the equation, 

-j— - t r t — Kx + Xot—X + Xott^—- 

orij dm orij on, orij 



(38) 



(39) 



In (39) for the sake of brevity we omit the spatial integrations. The variational derivative 
Sxo/Srii has the simple functional form, 



5x (r u r 2 ,uj) 
Srii 



[G (t u r 2 , to + Ei) + G (ri, r 2 , -u + e f )] <j>\ (ri)0i(r 2 ) 



(40) 



with G*o( r i) r 2? ^) being the Green function of iV non- interacting electrons moving in the single 
particle potential V# + v xc + V ext . We believe and some evidences can be found in jl| that 
the contribution coming from 5R/Srii is small. Nonetheless even having omitted this term, we 
are still to deal with rather involved eqs. (38,39). Therefore, it is of interest to illustrate the 
general consideration of the one-electron energies with a simple and important example, when 
only the exchange part A x of the total exchange-correlation functional is selected to be treated 
rigorously, using an approximation for A c = A xc — A x . Taking into account eq. (28) and the 
local density approximation, one gets, 



E c [p] = J p(r)e c ((r))dr. 



While for functional A x we have an exact expression [[05 



1 f \ dw 
A x [p] = E x [p] = -- / [xo(ri, r 2 , iw) + 27rp(r 1 )5(w)5(r 1 - r 2 )] f — dr 1 dr 2 . 

2 J \Ti — r 2 | 2n 

As it is seen from eq. (30), single-particle potential v xc takes from, 

v xc (r) = v x (r) + vJr), 



(41) 



(42) 



(43) 
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with the potentials being given by, 



SE X 
5p(r) 



v r \r 



5E C 
Sp(r)' 



(44) 



We remark that v x (r) can be calculated exactly 



proximations to the potential, see e.g. |27] 



, while there are quite suitable ap- 
Presenting single-particle potential v xc (r) in 
such a way, we simplify the calculations a lot, preserving the mentioned above condition 
v xc (r — > oo) — > — 1/r, which of crucial importance for calculations of £j related to the Kohn- 
Sham unoccupied states, see eq. (31). In the same way as eq. (38) was derived, one gets, 



< fa\v z 



_1 _5_ 

2 5n, 



Xo(ri, r 2 , m>) + 2tx p(r 1 )5(w)5(r 1 - r 2 ) 



ri - r 2 



2^ ' 



(45) 



After straightforward calculations one arrives at rather simple result for the single-particle 
spectra, that are to be compared with experimental results, 



6i = e~ < fa\v x \fa > - ^2n k 



|ri - r 2 | 



dvidv2- 



(46) 



The single-particle levels e^, given by eq. (46), resemble the ones that are obtained within 
the Hartree-Fock (HF) approximation. If the wave functions fa would be solutions of the HF 
equations and correlation potential v c (r) were be omitted, then the energies e« would exactly 
coincide with their eigenvalues. But this is not the case since fa are solutions of eq. (31), and 
the energies do not coincide with either HF eigenvalues or with the ones of eq. (31). We 
also anticipate that eq. (46) applying to solids will produce a quite reasonable results for the 
energy gap at various high-symmetry points in the Brillouin zone. 



CONCLUSIONS 



We have shown that it is possible to convert DFT into a self-consistent theory by constructing 
the exact equation for the exchange-correlation functional A xc . The proposed version of DFT 
has another desirable feature: a systematic way to construct successive approximations to A xc 
is becoming available. We have shown that the usual action functional A[p] can be used to 
calculate the response functions. It was also confirmed that the second functional derivative 
of the action functional with respect to the density defines the linear response function. The 
higher order functional derivatives define the higher order response functions. All they are not 
the causal ones, being symmetric under exchange of their arguments. We have also shown that 
on the same ground one can calculate the retarded (causal) and advanced response functions. 
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Thus, the paradox related to the symmetry properties of the effective interaction and the causal 
properties of the linear response function is resolved. The developed DFT allows for calculations 
of excitation spectra of any multielectron system. We have also related the eigenvalues of the 
single-particle Kohn-Sham equations to the real single-particle spectrum. 
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